#include "scene.h"
#include <stdio.h>
#include <math.h>


extern char **parse_str(const char *str, char c, unsigned short start, unsigned short end, unsigned short *nsubs);
extern double str2num(const char *str, short nsubs);

static char r_grid_gen(struct grid_config *gc, struct grid *g);
static char th_grid_gen(struct grid_config *gc, struct grid *g);
static char z_grid_gen(struct grid_config *gc, struct grid *g);

extern void grid_gen(struct grid_config *gc, struct grid *g, char *ecode)
 {
  ecode[0] = ecode[1] = ecode[2] = 'o';

  ecode[0] = r_grid_gen(gc, g);
  ecode[1] = th_grid_gen(gc, g);
  ecode[2] = z_grid_gen(gc, g);

  return;
 }


static char r_grid_gen(struct grid_config *gc, struct grid *g)
 {
  FILE *ptr;
  unsigned short i, cnt, isec, izone, nr, n;
  double *r;
  char fname[MAX_FNAME_LEN+1], *str, **s;

#ifdef READ_GRID_FILES
  sprintf(fname, "%02d.%s.dat", RUN_CNT-1, R_FILE);
  ptr = (FILE *)fopen(fname, "r");

/* Skip the 1st-line. */
  while(fgetc(ptr) != '\n');

/* Determining maximum line-length. */
  for(n = isec = 0; isec < Nsec; isec++)
    {
     for(i = 0; getc(ptr) != '\n'; i++);
     if(i > n) n = i;
    }

  str = (char *)malloc(sizeof(char)*(n+1));
  rewind(ptr);

/* Skip the 1st-line. */
  while(fgetc(ptr) != '\n');

/* Read the grid-points. */
  for(isec = 0; isec < Nsec; isec++)
    {
     for(i = 0; (str[i] = getc(ptr)) != '\n'; i++); str[i] = '\0';
     s = parse_str(str, ' ', 1, 2*NrG+g->nr[isec]+1, &n);

     if(n != 2*NrG+g->nr[isec]+1)
         {
          fprintf(g->lptr, "ERROR:: In 'r_grid_gen()': the r-grid file '%s' seems to be currupt.\n", fname);
          return 'r';
         }

     for(i = 0; i < n; i++)
        {
         g->r[isec][i] = str2num(s[i], 3);
         free(s[i]);
        }
     free(s);
    }

 free(str); fclose(ptr);

#else

# ifdef r_UNIF_GRID
  fprintf(g->lptr, "ERROR:: In 'r_grid_gen()': No rule defined for uniform grid generation in r-direction.\n");
  return 'r';
# else

/* Generating grids for the internal domains of all the sections. */
  g->r[0][NrG] = g->r[1][NrG] = gc[1].r[0]; g->r[2][NrG] = gc[1].r[5];

  for(isec = 0; isec < Nsec; isec++) for(cnt = NrG+1, izone = 0; izone < gc[isec].nzone_r; cnt += gc[isec].nr[izone++])
  for(i = 0; i < gc[isec].nr[izone]; i++)
   g->r[isec][cnt+i] = g->r[isec][NrG] + gc[isec].r[izone] + (i+1)*(gc[isec].dr[izone] + 0.5*i*gc[isec].r_cd[izone]);

/* For the ghost-cells. */
  for(isec = 0; isec < Nsec; isec++) for(i = 1; i <= NrG; i++)
     {
      g->r[isec][NrG+g->nr[isec]+i] = 2*g->r[isec][NrG+g->nr[isec]] - g->r[isec][NrG+g->nr[isec]-i];
      g->r[isec][NrG-i] = 2*g->r[isec][NrG] - g->r[isec][NrG+i];
     }
# endif

#endif

/* Printing the r-gid paramaters in a file. */
  sprintf(fname, "%02d.%s.dat", RUN_CNT, R_FILE);
  ptr = (FILE *)fopen(fname, "w");

/* Printing for backup if needed later in recovering a program. */
  for(isec = 0; isec < Nsec; isec++) fprintf(ptr, "%d ", g->nr[isec]);
  fprintf(ptr, "%d %d %d\n", g->nr_nozzle, g->nr_gap, NrG);

  for(isec = 0; isec < Nsec; isec++)
     {
      for(i = 0; i <= 2*NrG+g->nr[isec]; i++) fprintf(ptr, "%0.15E ", g->r[isec][i]);
      fprintf(ptr, "\n");
     }

  fprintf(ptr, "----------Don't edit above this line----------\n\n");

  for(isec = 0; isec < Nsec; isec++) fprintf(ptr, "# nr[%d] = %d\n", isec, g->nr[isec]);
  fprintf(ptr, "# nr_nozzle = %d\n# nr_gap = %d\n# NrG = %d\n", g->nr_nozzle, g->nr_gap, NrG);

  for(isec = 0; isec < Nsec; isec++)
       {
        r = g->r[isec]; nr = g->nr[isec];
        fprintf(ptr, "\n\n\t[ SECTION - %d ]\n# nr = %d\n", isec, nr);
        for(i = 0; i <= 2*NrG+nr; i++)
           {
            fprintf(ptr, "\n%d\t%E", i-2, r[i]);
            if(i < 2*NrG+nr)
               {
                if(i > 0) fprintf(ptr, "\t|------>%E", r[i+1]+r[i-1]-2*r[i]);
                fprintf(ptr, "\n\t\t|------>%E", r[i+1]-r[i]);
                if(r[i+1]-r[i] < 0)
                    {
                     fprintf(g->lptr, "ERROR:: In 'r_grid_gen()': r[%d] < r[%d].\n", i+1, i);
                     return 'r';
                    }
               }
           }
        fprintf(ptr, "\n");
       }
  fclose(ptr);

  return 'o';
 }


static char th_grid_gen(struct grid_config *gc, struct grid *g)
 {
  FILE *ptr;
  unsigned short cnt, i, n;
  char fname[MAX_FNAME_LEN+1], *str, **s;
  double *th;

 th = g->th;

#ifdef READ_GRID_FILES
  sprintf(fname, "%02d.%s.dat", RUN_CNT-1, TH_FILE);
  ptr = (FILE *)fopen(fname, "r");

/* Skip the 1st-line. */
  while(fgetc(ptr) != '\n');

/* Determining maximum line-length. */
  for(n = 0; getc(ptr) != '\n'; n++);

  str = (char *)malloc(sizeof(char)*(n+1));
  rewind(ptr);

/* Skip the 1st-line. */
  while(fgetc(ptr) != '\n');

/* Read the grid-points. */
  for(i = 0; (str[i] = getc(ptr)) != '\n'; i++); str[i] = '\0';
  s = parse_str(str, ' ', 1, 2*NthG+g->nth+1, &n);

  if(n != 2*NthG+g->nth+1)
         {
          fprintf(g->lptr, "ERROR:: In 'th_grid_gen()': the th-grid file '%s' seems to be currupt.\n", fname);
          return 't';
         }

  for(i = 0; i < n; i++)
        {
         g->th[i] = str2num(s[i], 3);
         free(s[i]);
        }

  free(s); free(str); fclose(ptr);

#else

# ifdef th_UNIF_GRID
  for(cnt = 0; cnt <= 2*NthG + g->nth; cnt++) g->th[cnt] = (cnt - NthG)*g->dth;
# else
  fprintf(g->lptr, "ERROR:: in 'th_grid_gen()': No definition for non-uniform grid in theta-direction.\n");
  return 't';
# endif

#endif

  sprintf(fname, "%02d.%s.dat", RUN_CNT, TH_FILE);
  ptr = (FILE *)fopen(fname, "w");

/* Printing for backup if needed later in recovering a program. */
  fprintf(ptr, "%d %d\n", g->nth, NthG);
  for(i = 0; i <= 2*NthG+g->nth; i++) fprintf(ptr, "%0.15E ", g->th[i]);

  fprintf(ptr, "\n----------Don't edit above this line----------\n");
  fprintf(ptr, "\n# nth = %d\n# NthG = %d\n\n", g->nth, NthG);
  for(i = 0; i <= 2*NthG+g->nth; i++)
           {
            fprintf(ptr, "\n%d\t%E", i-2, th[i]);
            if(i < 2*NthG+g->nth)
               {
                if(i > 0) fprintf(ptr, "\t|------>%E", th[i+1]+th[i-1]-2*th[i]);
                fprintf(ptr, "\n\t\t|------>%E", th[i+1]-th[i]);
                if(th[i+1]-th[i] < 0)
                    {
                     fprintf(g->lptr, "ERROR:: In 'th_grid_gen()': th[%d] < th[%d].\n", i+1, i);
                     return 't';
                    }
               }
           }
  fprintf(ptr, "\n");

  fclose(ptr);

  return 'o';
 }


static char z_grid_gen(struct grid_config *gc, struct grid *g)
 {
  FILE *ptr;
  unsigned short cnt, nz, n;
  short izone, i, isec;
  double *z;
  char fname[MAX_FNAME_LEN], *str, **s;

#ifdef READ_GRID_FILES
  sprintf(fname, "%02d.%s.dat", RUN_CNT-1, Z_FILE);
  ptr = (FILE *)fopen(fname, "r");

/* Skip the 1st-line. */
  while(fgetc(ptr) != '\n');

/* Determining maximum line-length. */
  for(n = isec = 0; isec < Nsec; isec++)
    {
     for(i = 0; getc(ptr) != '\n'; i++);
     if(i > n) n = i;
    }

  str = (char *)malloc(sizeof(char)*(n+1));
  rewind(ptr);

/* Skip the 1st-line. */
  while(fgetc(ptr) != '\n');

/* Read the grid-points. */
  for(isec = 0; isec < Nsec; isec++)
    {
     for(i = 0; (str[i] = getc(ptr)) != '\n'; i++); str[i] = '\0';
     s = parse_str(str, ' ', 1, 2*NzG+g->nz[isec]+1, &n);

     if(n != 2*NzG+g->nz[isec]+1)
         {
          fprintf(g->lptr, "ERROR:: In 'z_grid_gen()': the z-grid file '%s' seems to be currupt.\n", fname);
          return 'z';
         }

     for(i = 0; i < n; i++)
        {
         g->z[isec][i] = str2num(s[i], 3);
         free(s[i]);
        }
     free(s);
    }

 free(str); fclose(ptr);

#else

# ifdef z_UNIF_GRID
  fprintf(g->lptr, "ERROR:: In 'z_grid_gen()': No rule defined for uniform grid generation in 'z' direction.\n");
  return 'z';
# else
/* Grid-generation for the inner-cells of all the three sections of the domain. */
  g->z[0][NzG] = 0; g->z[1][NzG] = h0_CYL; g->z[2][NzG] = h0_CYL-h2_CYL;

  for(isec = 0; isec < Nsec; isec++) for(cnt = NzG+1, izone = 0; izone < gc[isec].nzone_z; cnt += gc[isec].nz[izone++])
  for(i = 0; i < gc[isec].nz[izone]; i++)
     g->z[isec][cnt+i] = g->z[isec][NzG] + gc[isec].h[izone] + (i+1)*(gc[isec].dh[izone] + 0.5*i*gc[isec].h_cd[izone]);

/* For the ghost-cells. */
  for(isec = 0; isec < Nsec; isec++) for(i = 1; i <= NzG; i++)
     {
      g->z[isec][NzG+g->nz[isec]+i] = 2*g->z[isec][NzG+g->nz[isec]] - g->z[isec][NzG+g->nz[isec]-i];
      g->z[isec][NzG-i] = 2*g->z[isec][NzG] - g->z[isec][NzG+i];
     }

# endif

#endif

/* Printing the z-gid paramaters in a file. */
  sprintf(fname, "%02d.%s.dat", RUN_CNT, Z_FILE);
  ptr = (FILE *)fopen(fname, "w");

/* Printing for backup if needed later in recovering a program. */
  for(isec = 0; isec < Nsec; isec++) fprintf(ptr, "%d ", g->nz[isec]);
  fprintf(ptr, "%d\n", NzG);

  for(isec = 0; isec < Nsec; isec++)
    {
     for(i = 0; i <= 2*NzG+g->nz[isec]; i++) fprintf(ptr, "%0.15E ", g->z[isec][i]);
     fprintf(ptr, "\n");
    }

  fprintf(ptr, "----------Don't edit above this line----------\n");
  for(isec = Nsec-1; isec >= 0; isec--) fprintf(ptr, "\n# nz[%d] = %d", isec, g->nz[isec]);
  fprintf(ptr, "\n# NzG = %d\n", NzG);

  for(isec = Nsec-1; isec >= 0; isec--)
       {
        z = g->z[isec]; nz = g->nz[isec];
        fprintf(ptr, "\n\n\t[ SECTION - %d ]\n# nz = %d\n", isec, nz);
        for(i = 2*NzG+nz; i >= 0; i--)
           {
            fprintf(ptr, "\n%d\t%E", i-2, z[i]);
            if(i > 0)
               {
                if(i < 2*NzG+nz) fprintf(ptr, "\t|------>%E", z[i+1]+z[i-1]-2*z[i]);
                fprintf(ptr, "\n\t\t|------>%E", z[i]-z[i-1]);
                if(z[i]-z[i-1] < 0)
                    {
                     fprintf(g->lptr, "ERROR:: In 'z_grid_gen()': z[%d] < z[%d].\n", i, i-1);
                     return 'z';
                    }
               }
           }
        fprintf(ptr, "\n");
       }
  fclose(ptr);

  return 'o';
 }
